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PROCEDURES FOR THE COMPUTATION OF UNSTEADY 
TRANSONIC FLOWS INCLUDING VISCOUS EFFECTS 


ABSTRACT 

Modifications of the code LTRAN2, developed by 
Ballhaus and Goorjian^, which account for viscous 
effects in the computation of planar unsteady 
transonic flows are presented. Two models are 
considered and their theoretical development and 
numerical implemenation is discussed. Computational 
examples employing both models are compared with 
inviscid solutions and with experimental data. A 
description for use of the modified code is provided. 
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I . INTRODUCTION 


Solutions of planar inviscid unsteady transonic flowfields about thin airfoils 
are commonly obtained by time integration of the differential equation 
governing the velocity potential function. This method is particularly 
attractive because it permits the treatment of nonlinear flow phenomena 
including irregular shock wave motion. Development of the LTRAN2 code by 
Ballhaus and Goorjian^ has made available an efficient time-implicit finite 
difference algorithm for obtaining solutions to the low-frequency small 
disturbance transonic potential equation. This code is currently employed for 
many unsteady applications and has evolved as a useful tool for performing 
aeroelastic calculations and flutter analysis. 

While the inviscid flowfields computed by LTRAN2 provide a reasonable physical 
description for a wide range of unsteady flow conditions, these results are 
not adequate when viscous effects are significant. Such situations can 
involve complex flow structures characterized by the following phenomena: 

1. Shock-boundary-layer interaction altering the shock strength and 
location, 

2. Effective camber modification due to differences in the 
boundary-layer displacement on the upper and lower airfoil surfaces, 

3. Displacement and camber effects of the near wake. 

A detailed description of these effects can be provided through use of more 
exact equations governing the flow. It is not the purpose here to provide a 
complete description, but rather to account for viscous effects only as they 
impact upon the unsteady airfoil surface pressure distribution. The objective 
is to modify the basic LTRAN2 code in such a fashion that the aforementioned 
phenomena may be treated without seriously degrading efficiency. Inherent to 
this approach are several simplifying assumptions which are thought to prove 
adequate for aeroelastic applications. 
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It is postulated that the viscous boundary layer is in instantaneous 
equilibrium with the unsteady inviscid flow. A simple order of magnitude 

analysis indicates that this assumption is valid if the reduced frequency of 

2 

the unsteady airfoil motion is small . As a consequence, a steady form of 
the viscous modelling may be applied in a quasi-steady fashion. 

The most complex aspect of viscous transonic flows is that of the 

shock-boundary-layer interaction. A number of calculations for both steady 

2 3 

and unsteady flows ’ have indicated that this effect may be modelled by a 
simple computational artifice. This is accomplished by augmenting the surface 
geometry with a wedge-nosed ramp that is inserted at the base of the shock in 
an inviscid calculation. 

The ramp weakens the strength of the shock, usually displacing it upstream 
from its inviscid location, thereby simulating the effects of boundary-layer 
displacement. The ramp angle is chosen such that an empirically defined 
post-shock pressure is recovered. 

For many flows of interest, use of the ramp model alone will suffice to 
predict unsteady pressure distributions. When aft viscous effects are 
significant, the ramp is used in conjunction with an integral boundary-layer 
calculation over the rearward portion of the airfoil and the downstream wake. 
For this purpose "lag entrainment equations" ’ are employed, which have 
been found useful for the prediction of turbulent shear layers. In addition, 
a coupling of these equations with an inviscid algorithm has previously been 
used to compute steady transonic flowfields about airfoils, both with and 
without separated regions^’^. 

The sections that follow will first briefly review the theoretical development 
of the governing equations used to modify the LTRAN2 code. Numerical 
implementation and details of the modifications will be described. Finally, 
several computational examples will be presented and comparisons of the 
modified results will be made with inviscid solutions and with experiment. 
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II. THEORETICAL DEVELOPMENT 


The basic unsteady small disturbance perturbation potential equation treated 
by LTRAN2 may be written as 


,2kM^ , 

( ») 


a 


Z7T 


xt 



- {y + Dm"*^ 


<’x]*xx 



(1) 


where k = (oc/LI^ is the reduced frequency, y the ratio of specific heats, 
M ^ the freestream Mach number and subscripts denote partial 
differentiation. Here x and y are the streamwise and normal Cartesian 

-1/3 

coordinates normalized by c and c 5 respectively, where c is the airfoil 
chord and 5 the thickness ratio. The time, t, has been normalized 

by<o"^, u) being the circular oscillation frequency. The perturbation 

velocity potential function, <>, has been normalized by cU 5 ^^^ where 
is the freestream velocity and mis commonly taken to be 1. If the airfoil 

surface is defined by y" = “(x,t) the flow tangency condition may 

be expressed as 

()) " = f " on y = 0 “ for 0 < x £ 1 . (2) 

y ^ 


Across the trailing vortex wake the contact conditions 



(continuity of slope) (3a) 


on y = 0, X > 1 

A 


[ 



= 0 


(continuity of pressure) (3b) 
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are applied where the brackets denote the jump in the enclosed quantity from 
above to below the vortex sheet. At the outer boundaries the following 
conditions are applied: 

()) = 0 far upstream (4)- ' ‘ 

<)) = 0 far laterally (5) 

^ n 

= 0 far downstream. (6) 

X ^ 

h 

Finally, the initial condition (]> (x, y, 0) = I(x, y) is specified. 

A. Viscous Ramp 

The basis for development of the viscous ramp model is the observation 

in many steady experimental measurements that the post-shock pressure 

for turbulent flow over an airfoil corresponds approximately to that of 

the oblique shock produced by flow over a ramp with a wedge slope equal 

to the "detachment" angle . This observation has proven extremely 

useful as a simple computational artifice which accounts for the gross 

dominant effects of the shock-boundary-layer interaction. Computational 

results generated by inviscid codes which have been modified to 

incorporate the ramp model have been shown to correlate well with 

experimental data in steady, and more recently unsteady calculations 

2 8 

where the quasi-steady assumption has been invoked ’ . 

The geometry of the viscous ramp model is depicted schematically in 

Figure 1. It consists of a short precursor over which the surface slope 

varies from zero to the given wedge angle followed by the main ramp body ^ 

along which the slope varies quadratically. At the ramp leading edge 

the ramp height and slope are continuous, and at the downstream end the 

slope and curvature vanish. The ramp is positioned with respect to the 

local instantaneous sonic point location and affinely scaled with the 

shock strength as determined by conditions immediately ahead of the 

3 . • • ' , 

shock. As was done in the steady case , the ramp leading edge is " 
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offset a distance ahead of the sonic location in order to more 
fully influence the numerical shock profile. The leading edge is 
preceded by the precursor of length Xp which has been used in previous 
unsteady calculations^’® to moderate the passage of the sharp ramp 
leading edge across computational mesh points, thereby precluding 
spurious numerical instabilities. The main body of the ramp has a 
length of Xj^ where the parameters Xp, x^, and Xj^ may be selected 
for each specific application. 

Augmenting the surface geometry by the ramp model results in the 
following modification to the tangency condition, Eq. 2 : 

= f * + on y = 0 “ for 0 < x < 1 . (7) 

y X - - 

If X is the sonic point location and e the wedge angle then 

+ ^ 

fj^ is defined as follows: 


fp- 

= 0 for X < X - x„ - x„ 

(8a) 

®x 

sop 


p X ""X ^*X ’ T 

(8b) 

- Xq - Xp < X < X^- Xp 


X*X X 

fR* ' * ‘j-’ [l - • - Xo i i *5- Xo* Xr (80 

X K 

fjT = 0 for X > x^ - x„ + Xn . (8d) 

S 0 K 
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where 


The instaneous wedge angle 0, corresponding to that for the maximum 
possible turning for an attached oblique shock is found from the 
following analysis. Consider a wedge of angle 9 with an oblique shock 
located at its vertex. Within the approximations of small disturbance 
theory it is easily shown from the tangency condition that 


5 (i'y2 ■ 


= tan 9 9 


( 9 ) 


where the subscripts 1 and 2 denote the upstream and downstream states 
respectively. We now write Eq. 1 in conservative form as 






b (j) + c 


Vx 


u ) 

y y 


( 10 ) 


For convenience the following definitions have been made: 


k rf 


2/3 


b = 


1 - 
~27T 


(y + DM' 


c = - 


( 11 ) 


If i^. and i^ are unit vectors in each of the coordinate 
directions, the following mass flux vector, q, may be formed: 


q = (a + (a 4.^ + b 4.^ + c 4>j^^)i\ + ^‘*‘y^\ . 

Given the unit vector normal to the shock, IT, the shock jump conditions 
are then defined by 

[■q ]. T =0 (13) 

[ ]={ U - 
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Let us now assume that = x - g (y, t) = 0 describes the shock 
geometry. It follows that n is equal to, . 

i^h^i 


Eq. 13 then becomes 

- [a +^]gt + [a ft + + c [♦y] = 0 . (14) 

In a similar manner one may consider the jump conditions for the 
following relationships 



U),=o 

(15a) 

y X 

x y 

'♦x't- 

4,)x-0 

(15b) 


and thereby eliminate the shock geometry, i.e., g^ and g . Upon so 

^ y 

doing, one arrives at the "Rankine-Hugoniot" condition corresponding to 
Eq. 1: 


^a(*t2 ■ " *xl' ^ ^*x2 ’ 


(16) 




. ^^2 

For, maximum turning.-^ = 0, a condition which is obtained from 


34. 


x2 
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Eq. 16 as 


2b ♦xl 
*x2 " ■ 3c ■ 3 


(17) 


using 4>t2 " *tl ° • 

This relationship is now used in conjunction with Eq. 16 and Eq. 9 to 
find that 


9 

5 


m 


1-M" 


4 

^1 - -jrr- 

5 

(y+1)M 

— oo — ' 

3 


(18) 


Here the plus sign is taken for the upper surface and the minus sign for 
the lower surface. 

B. Lag Entrainment Equations 

The form of the so-called "lag entrainment equations" which will be 

considered here are those due to Green . They are predicated upon the 
boundary-layer assumption that the normal extent of the viscous region 

is small when compared with airfoil or wake thickness, which necessarily 

applies to flows at high Reynolds numbers. By integrating the governing 

partial differential equations in the normal direction and suitably 

modelling the requisite relationships among the dependent variables, a 

system of three first order ordinary differential equations is 

obtained. Two of the equations result directly from continuity and 

momentum. The third evolves from the Bradshaw, Ferriss and Atwell^^ 
turbulent kinetic energy equation but is formulated in terms of the 

entrainment concept originally proposed by Head^^. This yields a 
streamwise rate equation governing the degree to which the outer 

inviscid flow merges with the boundary layer. 


11 



The choice of an integral set of equations as opposed to a differential 

set is made for two reasons. First, integral techniques generally 

require considerably less computing time than differential methods. 

Second, differential methods have not as yet proven to be substantially 

superior for the computation of flows of practical interest. Since 

their original formulation the lag entrainment equations have proven to 

4 9 

be a useful computational tool for performing compressible ’ and 
incompressible boundary-layer and wake calculations, involving both 

c 7 1 p c 

separated ’ ’ and unseparated flow regions. 


For completeness. Green's lag entrainment integral equations are briefly 

summarized here. If a more thorough understanding of their development 

is required, the reader is referred to any of the original 
4 5 9 

references ’ ’ . It is useful to first note some prel imi nary 
definitions: 

* ” pU 

displacement thickness « = / (1 - pgU )d n (19) 

0 ® 

CO pU _U 

momentum defect thickness e = / (1 -1)^) d n (20) 


shape factor 




( 21 ) 


shape parameter 


= J J ^(1 --L) dn 
® 0 ^e ^e 


( 22 ) 


1 d ^ bl 

entrainment coefficient = — tT" f pUdn 

E PgUg d^ ^0 


(23) 
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In the preceding equations 5 and n are streamwise and normal physical 
coordinates respectively, p, the density and U the streamwise velocity 
component in the boundary layer. The subscript e refers to the 
boundary-layer edge denoted by n = Within the confines of small 

disturbance theory and the boundary-layer limit, the edge conditions may 
be taken as those corresponding to the inviscid solution on the airfoil 
surface or the wake centerline. 


The primary dependent variables are taken as 0 , H, and C^. Given 
their value at any streamwise station, 5 , their distributions may be 
predicted by the following first order set of ordinary differential 
equations and subsequent parameters: 


da ^f 2 a *^^e 


(24) 




(25) 


dC 


E=(L)[ {^) he )l/2.,c 

9 ^ H+H^ L ^ EQO ^ J Ug 


dU. 


d 5 'EQ 


- [. .0 

U, d; [- 


075 M 


Y -1 2 

y (I ^~T~ rW J- 
c 

® 1 + 0.1 


(26) 


The following relationships are now used to complete the description of 
the governing equations: 

Fc = [1 + (^) , (27) 
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locltn:' 


Fp = 1 + 0.056 Mg^ , 


(28) 


Re„ 


p U 9 
e e , 


(29) 


' = 1 - r 0.01013 

L logio(FRR,J 


- 0.00075 


- 1.02 


]- 


(30) 


. 9 -.1/2 

= TT { 1 - 6.55 [ ( — (1 + 0.04 Mg*^)] } , 


(31) 


H “1 

0.9 ( 0.4) - 0.05 for wall boundary layers 

TO u 

C. = J % 




0 for wake shear layers, 


(32) 


H = CH ■+ 1) [ 1 + (^ )rMg^ ] -1 , 


(33) 


H. = 3.15 + {^-^) - 0.01 (H - 1)^ , 
1 F-1 


(34) 


dH 




(H - 1)^ 


1.72 + 0.02 (F+1)' 




(35) 


0.8 C 


F = 


0.02 + 

0.01 + Cr 


fo 


(36) 


14 , 



( 37 ) 


C = (1 + 0.1 M^^)(0.024 C. + 1.2 + 0.32 ) . 

T e h L TO 




1 for wall boundary layers where secondary influence of 
X = i streamwise wall curvature have been neglected (38) 

1/2 for wake shear layers. 


(?r- ^)eo = ) (1 + 0.04 (39) 

H L ^ 6.432 TT ® -I » 


‘U, dc 


«=e'eqo“"i 


e ‘^^e 1 

^^eq] ’ 


(40) 


"=,'eQ0 = ‘1 * ( “•‘>2“ <'=E>EQ0 * [‘Ce’eQoJ 


+ 0.32 C 


fo 


(41) 


C = (C^)eqo(I + 0-1 - 0*32 , 


(42) 


(Q ) = ( ^ )^^^ - 0 01 

eq 1.2 + 0.0001 ' 


(43) 


dU 

(? 1) 

dc ^EQ 


[,g,>W][ V - 'vj ■ 


(44) 


Using the small disturbance approximation, edge qualities may be 
written in terms of the inviscid solution as follows: 


^e . ^ 2/3 
U " ^ 5 ‘*’x ’ 


(45) 
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(46) 


= H, ( 1 + [l * (1^) M„2] «2/3 } , 


P « X 


(47) 


t 

The pressure gradient parameter may be approximated as 


1 ^^e 1 2/3 ^ 

Ug (t c ^ Sx 


(48) 


and the local Reynolds number based on momentum defect thickness by 


Re, 






(49) 


where R_ = ^ o.^ «,c 


is the freestream Reynolds number based upon chord. 


The viscosity ratio ^ is now eliminated from Eq. (49) by noting that 

i 

— = 1 - ( Y - 1) 5^/3 M ^ A (50) 

^ X 

00 
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cr 



and then employing the Sutherland viscosity law^^ in the form 


(— ) = (-^) 
U_ T_ 



( 51 ) 


where So is a constant dependent upon the nature of the specific gas 

(So = 110°K, for air). It is noted in Eq.'s 26, 27, and 33 that r is 
the recovery factor which is commonly chosen as a function of the 

turbulent Prandtl number Pr^ : 


r = (Pr^) 


1/3 


(52) 


With these specifications the lag entrainment equations are fully 
defined. Noting that the shape parameter H and the entrainment 
coefficient are already nondimensional , the equations are rewritten 
here in a form compatible with coupling to the inviscid solution; 


— (-) = 
dx ^c' 


dx dH^ 


- (H + 2 - m/) ,9.) _ 

(Cj - Hj [^ - (H + 1) 5^^^ (|) 




(53) 


(54) 
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- 6 


2/3 /9 


C XX 


(1 


1 + 0.075 M. 




2^. 


1 + 0.1 M 


Me^) 

T 


(55) 


A — 

Given initial values of (— ), H, and at the streamwise station x = x^. , 
Eq.'s 53 through 55 may be integrated in the downstream direction for 
all X > x^. to obtain the distribution of viscous parameters along the 
airfoil surface and trailing wake. 

These are then used to modify the inviscid solution through use of the 
displacement thickness, s*. An explicit dependence of the slope of 

ds * 

the displacement thickness, upon the viscous parameters provided by 
the lag entrainment equations may be obtained by the following 
analysis. From Eq. 21 we write 


ds* 

dc 



dH 

^ • 


Differentiating Eq. 33 results in the expression 


(56) 


(k 


= [ 1 + (^^)r ^ + (y-1) r Me (H + 1) ^ . (57) 


Now, using the isentropic relationship 


._13 


Me 


[l * (lii) ] 


- 1/2 


(58) 


where a^ is the speed of sound at the stagnation condition, it may be 
shown that 


M 

k e 


dU, 




(59) 
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Upon substitution of Eq. 59 into Eq. 57 to obtain 


dH _ f, . M 21 dH 

^ - [1 + (T-)r J ^ 

t (r-l)r h/ 1 + (1'^) 


( 60 ) 


and subsequent use of Eq. 60 in conjunction with Eq. 56 results in 


ds* _ „d9 

ar ■ w 


+ 9 




+ (•,-!) (F-H)[l + (1^) m/ ] (61) 


Finally, the derivatives on the right-hand side of Eq. 61 are eliminated 
by employing Eq.'s 24 and 25 to arrive at 


W ^ + Ll + (V> >• «e J I'^E - 

+ (v-l)r (H+1)] - H (H+2 - M^^)l (^|p) • (62) 
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By nondimensional izing Eq. 62 and applying the small disturbance 
assumption, we can transform this relationship to a form compatible with 
the inviscid equation as 


d 

■31 



Fi .F2 


XX 


where for convenience we have defined 


(63) 


= { ^ + [1 + (^)r ](C^ - -^^1/ 5'\ (64a) 

^2 " ^ ^ ^ ^e^][^l (H+l)](64b) 

- H (H+2-M ^)} (-) . 

e c 


5 

The form of Eq. 63 was first inferred by East 
interesting properties. For unseparated flows, 


et. a1 
F 2 4^ 


. who noted several 


XX 


<< 


d 5 

so that Eq. 63 represents a stiff equation for ^ (— ) . 

Thus, the displacement is relatively insensitive to the pressure 
gradient, <{ 1 ^^. On the other hand, for fully separated flows 

X X 



>> 



In this case Eq. 63 may be rewritten as 




• which is a stiff equation for Eq. 63 is now used to provide 
the coupling between the inviscid and viscous solutions. The original 


i 
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surface boundary condition, Eq. 2, is modified by the viscous 
displacement thickness to yield 

♦ = f * + (|^)^ on y = 0- for 0 < X < 1. (65) 

y A 0 c X — — 

Using Eq, 63 this is written as 

<t»y" = + (Fj^ + p 2 on y = 0* for 0 £ X _< 1 (66) 

in order to exhibit the dominant explicit dependence of the displacement 
upon the velocity potential. 

Downstream of the trailing edge, the viscous wake generates an effective 
displacement afterbody. Due to its presence, a discontinuity in the 
slope of the potential is now permitted such that the wake condition Eq. 
3a must be replaced by 


[■fry] ^ [ ^1 *^2 *xx] on y = 0 for X >1 . 


(67) 
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III. NUMERICAL METHOD 


A. Viscous Ramp 

Given the velocity potential distribution 4 >(x, y, t*^) over the 
computational domain at time level t*^, the viscous ramp model is 
implemented in the following manner in order to extend the solution to 
time level = t*^ + At. First, the sonic point, x^, defined as 
that location where the surface pressure attains the critical value 
(Cp*), is located based upon ())(x, y, t"^). A simple linear 
interpolation between mesh points is employed for this purpose. Next, 
the condition ^ is evaluated where the upstream location is taken 
to be that position which lies two mesh points upstream of the first 
subsonic point. The first subsonic point is that mesh point immediately 
downstream of x^. With now known, the value of the ramp angle 

is calculated from Eq. 18. Given this value and the sonic point 
location, the ramp geometry is computed according to Eq. 8. Finally, 
the modified boundary condition, Eq. 7, may be applied and the solution 
advanced to the new time level 


The ramp geometric parameters Xp, x^, and Xj^ can be adjusted for 
specific applications to provide desired computational results. This is 
often a matter of trial and error and it is suggested that comparisons 
with steady experimental pressure distributions, if available, be used 
for this purpose. The modified LTRAN2 code will accept these variables 
as input parameters. In addition, the following program default values 
have been proven adequate for a wide range of both steady and unsteady 
calculations: 


0 . 02 , 

0 . 02 , 

0 . 10 . 
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B. Lag Entrainment Equations 


The coupling of the lag entrainment equations with an inviscid solution 
algorithm can be an especially sensitive procedure. This is 
particularly true in the transonic regime where it is in general 
necessary to integrate the integral equations through the numerical 
shock profile. The transition from supersonic to subsonic flow regions 
can often produce numerical oscillations, instabilities, or other 
anomalous behavior. In addition, the structure of the shock profile is 
typically determined by the numerical techniques of the shock capturing 
method rather than a strict modelling of the physical processes 
involved. For example, the use of backward and centered streamwise 
differencing in the supersonic and subsonic regions respectively 
generally produces a shock profile with a streamwise extent of two mesh 
intervals. Thus the nature of the solution in this region is very much 
dependent upon the mesh distribution unless it is exceedingly fine. 
Although these deficiencies are not serious enough to degrade the 
accuracy of the results for aeroelastic applications, any rigorous 
coupling of the viscous and inviscid solutions through the shock profile 
is necessarily superfluous. 

For these reasons, use of the viscous ramp model is maintained in 
conjunction with the lag entrainment equations. As was done previously, 
it is assumed that the solution (|){x, y, t*^) at time level t^ is 
known over the entire computational domain. It is consistent with the 
quasi-steady assumption, then, to use this solution to integrate the lag 
entrainment equations and extend the potential distribution to the next 
time level t"^^ = t"^ + At. An initial upstream position, x^. , is 
selected as a location ahead of the shock. It is assumed that upstream 
of x^. the effect of the boundary layer will be negligible. 



Initial values for the quantities (— ), H, and Cr at station x. are 

V W I 

obtained from empirical correlations for the incompressible turbulent 
boundary layer on a flat plate^^. In particular, we assume 

~ = 0.04625 X Re . (68) 

C 00 


Then using the relationship 


L. - 
■ Si 


(69) 


and Eq.*s 19 through 23 we can find that 


c 


(- 1 ) > 


(70) 


H = 



(71) 


d 5 * 

C = 7 — (M . (72) 

E dx c 

With the initial conditions specified, the parameters So, Re^, , 

and Pr. defined, and the distribution i))(x, y, t^) given, Eqs. 53 
through 55 may be integrated from the station x = x^- + ax to the 

downstream boundary. For this purpose a second order accurate 

Runge-Kutta predictor-corrector numerical integration scheme was 

employed. While other more accurate techniques could easily have been 

selected, the truncation error incurred in the chosen method is of the 

same order as that of the inviscid algorithm. The evaluation of the 

expressions in Eqs. 45 thorugh 49 was performed using a central 
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difference to calculate i{> and either a central or backward difference 

A 

to compute based upon the local flow characteristics (i.e., a 
central difference was used in subsonic regions and a backward 
difference was used in supersonic regions). This treatment is not only 
consistent with the inviscid algorithm, but also duplicates the 
technique employed in previous steady computations®’^. Once the 
integration has been completed the functions and F 2 defined in 
Eq. 64 may be evaluated along the airfoil surface and entire trailing 
viscous wake. These are then employed to couple the viscous soluion to 
the inviscid solution through the use of Eqs. 66 and 67. 

Although integration of the lag entrainment equations is initiated 
upstream of the shock, coupling with the inviscid flow does not occur 
until the first mesh point downstream of the shock for which the flow is 
locally subsonic, defined as This ensures that the turbulent 

variables have properly adjusted to the shock-boundary-layer interaction 
when the coupling procedure is invoked. Upstream of the viscous 

ramp technique is implemented in the fashion described in the previous 
section. Beginning at Eq. 66 is used as the modified boundary 

condition for the inviscid solution. This treatment can potentially 
produce anomalous behavior due to inconsistencies between the respective 
methods, but none have been observed and the procedure appears to be 
reasonably stable numerically. 

Equation 66 is approximated by evaluating 4 >y as a first order 
one-sided difference and <t>xx ^ second-order central difference, as 
the condition is only applied in subsonic flow regions. If it is 

n+1 

assumed that the solution at the new time level t" has been 
generated for all x <x j, then Eq. 66 for the upper surface may be 
expressed at x = Xj as follows: 



n+1 n+1 \ 

’I, JW+1 " ’I, JW \ ^ ^ c 

W- <ow 



( 73 ) 


where yjy = 0’*’. It is noted that at the station just downstream of 
the station under consideration, a value for ^ at the previous 

time level must be taken in order that the inviscid algorithm remain 

intact. Equation 73 is now rearranged such that all known quantities 

appear on the right-hand side: 


(y Y ^ 

'^JW+l ■ '^JW 


n+1 

I, JW+1 



1 

'jw+l ’jw 


n+1 


I, JW 



*1+1, JW 



(74) 


When used in conjunction with the difference expressions for the 
flowfield differential equation, this relationship provides an implicit 
expression simultaneously coupl ing all points above the airfoil 

(y ^ 0^) at the station x,. This treatment is considered to be 
significant in suppressing numerical instabilities which would otherwise 

occur, particularly in regions of large pressure gradient such as those 

at the trailing edge. An expression similar to Eq. 74 may easily be 

obtained for the lower surface. 
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Downstream of the trailing edge, the effects of the viscous wake are 
manifested in the jump condition, Eq. 67. This impacts upon the 
inviscid algorithm through the term in the governing differential 
equation when evaluated along y = 0 for x >1. It follows from Eq. 36 
that once the jump in potential at the trailing edge is found at time 
level t”'*’^, it is uniquely determined for all subsequent downstream 
stations. Thus (t.^^ may be evaluated along either y = 0^ or y = 0" 
and the value along the opposite side of the wake determined from Eq. 36. 

We now consider (4«yy)j where yj^ = 0" and x > 1 and for 
convenience define 

_ n+1 _ . /V n'*' . /V n~ 

Fj - 4* 'Xj, 0 , t ) ~ <t> (Xj, 0 , t ) , 

3j = ♦y (Xj, O'", t"^^) - ♦y {Xj, o', t"^^) . 

Using Eq. 67, we can express Eq. 76 as 

■ *" 11 ^ ^'"21 *xxl ■ ^"21 *xxl^ * 

Because the wake is entirely subsonic, a central difference is used to 
approximate the derivatives in Eq. 77 such that 


(75) 

(76) 


(77) 
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Values of the potential on the upper side of the wake are now eliminated 
through use of Eq. 75. In addition, we note that the potential at the 
upstream station Xj_j^ at the new time level are known, those at the 
current station Xj are to be solved for, and those at the downstream 
station Xj^.j must be taken at the previous time level t*^. Thus we 
may write 
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(79) 


A second order approximati 


on for j„ 


may now be written as 


/l, JW+1 JWx ,’I, JW JW-1^ 

( — y ^ ) -( — ^ ) 

,, , JW+1 JW JW JW-1 ,on^ 

‘‘"yy'l.JW" ^ ^ ^ 


'^JW+1 ■ ^JW- 




where (|>j is a modified value of 4>(Xj 
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Because 4. and ^ are discontinuous on y = 0, the value of 

-s» ^ 

((>1 must be computed in terms of 4>j jw+i» ^i» sndBj. 

From a Taylor series expansion about y = 07 we write 


*I,JW+1 ^ ^I, JW ^^JW+1 “ *YI ^ L 

Similarly, expanding about y = 0+, 

^I,JW+1 ^^JW+1 " *YI ° t^^JW+1 “ 

Upon combining Eqs. 81 and 82 along with Eqs. 75 and 76, the following 
expression is obtained: 


*1, jw+1 " '^i, jw+r ■ ^^jw+i ■ • 


(83) 


Finally, Eq. 83 is substituted into Eq. 80, the resultant expression 

evaluated at time level n + 1 and Bt eliminated by use of 
Eq. 79 to arrive at the expression 
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Equation 84 now provides the implicit condition coupling the inviscid 
algorithm with the viscous downstream wake. 
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IV. COMPUTED RESULTS 


In order to validate the modifications to the LTRAN2 code and to establish the 
acceptability of generated solutions, several computational examples of 
practical interest were considered. Two airfoil sections were selected as 
being representative of those for which viscous effects may be significant, 
namely the RAE 2822 and the NLR 7301. Section geometries are shown in 
Figure 2. The RAE 2822 section has a maximum thickness to chord ratio of 
0.1210 and is characterized by appreciable camber, but a thin nose region. On 
the other hand, the NLR 7301 section has a very thick blunt leading edge, a 
large aft cove, and a thickness ratio of 0.1683. These sections, particularly 
the NLR 7301, are considered to be a severe test of the limitations inherent 
in the computational method. This selection is based upon the existence of 
experimental data available for comparative purposes^^’^®. 

All of the results presented here were generated on a nonuniform Cartesian 
grid consisting of 113 points in the x-di recti on and 97 points in the 
y-direction, symmetric about y = 0. Minimum grid spacings were aXj^^^ = 

0.005 at the leading edge, 0.01 at the surface, and the computational 

domain was defined by 


-200 ^ X ^ 200, 

-397.8 < y < 397.8. 


The airfoil surface was defined by 48 points in the x-di recti on. For each of 
the cases considered, results were generated for an inviscid solution, the 
viscous ramp solution, and the lag entrainment solution. Comparisons between 
these results, and with experiment when available, were then made. 

A. The RAE 2822 Airfoil 

Solutions for the steady flow at M = 0.73, a = 3.19°, 

6 ™ 

and Re = 6.5 x 10 over the RAE 2822 section were obtained by 

no 

integrating in time. The inviscid solution was evolved from an initial 
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undistributed flow (i.e., <i) =0). This result was then used as the 
initial flow to obtain the viscous ramp solution, which in turn provided 
the initial condition for the lag entrainment result. While this 
process may not in general be optimal, it was employed here because all 
three solutions were required for comparison, and has been used to 
generate all subsequent steady results. Implementation of either the 
viscous ramp or lag entrainment computation (which is used in 
conjunction with the ramp) usually requires a local minimum of the 
pressure occurring over the airfoil surface in order to initiate the 
computation. This provides a starting location upstream of the trailing 
edge for the viscous effects to impact upon the solution as was noted in 

3 

earlier formulations for steady solutions obtained by relaxation . 

Surface pressure distributions for this case appear in Figure 3 and are 

15 

compared with the experimental results of Cook et al. It is seen 
that the inviscid shock location lies downstream of the experimental 
position. The inclusion of viscous effects displaces the shock 
forward. Location of the shock appears to be predicted quite well by 
use of the wedge alone, but not by the lag entrainment solution. The 
reason for this is that this case is particularly sensitive to small 
perturbations in the flow conditions. Other computations^^, for 
example, have slightly varied the Mach number and angle of attack to 
obtain better agreement with the experiment. Just upstream and 
downstream of the shock it is noted that the lag entrainment solution 
produces a pressure distribution that has no appreciable gradient. This 
evidences the fact that the shock position has little perf erred location 
and will be quite sensitive to small perturbations in the flow 
conditions. The lag entrainment solution is also seen to produce better 
agreement with experiment upstream of the shock, on the lower surface, 
and in the trailing edge region. In this result the viscous 
displacement has properly accounted for the effective decambering of the 
airfoil. Because this effect is absent in the wedge solution, its close 
agreement with the experimenal shock location may be fortuitous. 
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Distributions of the upper surface displacement thickness and momentum 
defect thickness corresponding to this case appear in Figure 4. The 
comparison of the computed solutions with experiment is qualitatively 
correct. Downstream of the shock, the levels of the computation lie 
somewhat below the experiment. This is probably because the shock 
location is poorly predicted. If the shock were located further aft, 
the rise of 9 and s* through the shock would be expected to be greater 
due to a larger local Reynolds number in the wake region. In addition, 
the behavior may be due in part to the fact that the inviscid and 
viscous solutions are not coupled through the shock profile. Downstream 
of the trailing edge 0 and 6* can be seen approaching a constant value 
which is consistent for free shear layers without pressure gradient. 

This is indicated, for example, in Eq. 53 with = 0. 

In order to exemplify the sensitivity of this computation to the flow 
conditions, a second solution was generated for = 0.75. Pressure 
distributions for this case are shown in Figure 5. Upon comparing this 
result with those in Figure 3 for M ^ = 0.73, it can be seen that the 
slight increase in Mach number has produced a significant effect on the 
numerical solutions. It is observed in the inviscid computation that 
the shock has moved to the trailing edge and that in the wedge solution, 
the shock location is far aft of the experimental position. The lag 
entrainment result, on the other hand, compares quite well with the data 
and appears to duplicate a previous steady calculation^. Once again, 
the reason for this is because the lag entrainment solution accounts for 
the aft displacement effects, while the other results do not. This 
effect is also noticeable on the lower surface and in the regions both 
ahead of the shock and at the trailing edge. The anomalous behavior of 
the wedge result just upstream and downstream of the shock probably is 
produced by the positioning of the wedge within the shock profile. This 
behavior could possibly be eliminated by alternate choices for the wedge 
parameters (Xp, x^, X[^) , but no attempt was made to do that here. 
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Corresponding upper surface distributions of the boundary layer 
quantities are presented in Figure 6. The lack of agreement with the 
data occurs because the flow was found to separate experimentally, but 
remained attached for the computed solution. Much higher levels in 9 
and 5 * naturally occur downstream of separation. This condition is once 
again probably quite sensitive to the flow conditions. It is 
interesting to note that the numerical results of Reference 7 
overpredicted the level of the 5 * distribution, while those presented 
here have underpredicted it. The fact that the pressure distributions 
of both results compare well with experiment is due to the fact that the 
coupling between the inviscid and viscous solutions is provided by 

d5* which is insensitive to ik as was noted earlier. 

W 

As an example of an unsteady calculation, this airfoil was forced to 
oscillate sinusoidally in rotation about the mid-chord at M = 0.73 

00 

with a reduced frequency k = 0.2 and an amplitude of 1°, The solution 
was initiated from the steady-state results in Figure 3 and the 

oscillation was allowed to proceed for three cycles with At selected 

such that 0.5° of rotation occurred in one time step. Figure 7 shows 
the instantaneous surface pressure distributions for the last cycle of 

oscillation where the time is referenced to the beginning of this cycle 

and given in degrees. The inviscid result for this case could not be 

computed with the same time step, and thus does not appear in the 

comparison. 

As was noted in the steady pressure distributions, the lag entrainment 
solution produces a weaker shock with more forward displacement than 
occurs for the wedge alone. Also, the effects of decarabering on the 
lower surface and in the trailing edge region are noticeable. At 
t = 360 the wedge is seen to generate a post shock re-expansion. Use of 
the viscous ramp can cause this behavior which may be altered by 
variation of the ramp parameters. 
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Instantaneous upper surface distributons of the displacement and 
momentum defect thickness for this case appear in Figure 8. They 
exhibit the same general properties of the steady-state results 
(Figure 4). In Figure 9 the time histories of the unsteady lift and 
moment coefficients are presented. Here the moment coefficient is taken 
about the pitch point x = 0.5. The unsteady angle of attack is provided 
for comparison. Phase differences between the lift and moment and 
oscillatory motion can be clearly seen. 

B. The MLR 7301 Airfoil 

Flow about the NLR 7301 section provides a severe test for the 
computational method described in preceding sections. Not only is the 
airfoil exceptionally thick, but the leading edge is extremely blunt. 
Thus, the assumptions inherent in small disturbance theory may be' in 
violation. Nevertheless, this case was selected because of available 
existing unsteady experimental data. In order to render the geometry 
more tractable for a small disturbance calculation, the original airfoil 
slopes were slightly modified. This modification is discussed in 
Appendix A and has been used for all the computations presented here 
with the NLR 7301 airfoil. 

Steady state solutions were generated for M = 0.75, a = 0.37°, and 
7 

Re = 1.148 X 10 in the manner as was described for the RAE 2822 

00 

airfoil. These conditions duplicate those of the experiments of Davis 
and Malcolm^^. Surface pressure distributions for the respective 
results and the experimental data appear in Figure 10. Both the 
inviscid calculations and the wedge solution are seen to agree favorably 
with experiment upstream of the shock. Location of the shock is poorly 
predicted by the inviscid result. Without the aforementioned slope 
modification, the shock was found to lie at the trailing edge. 
Improvement occurs when the wedge is included. On the lower surface 
both solutions are seriously deficient. The lag entrainment result 
appears to predict the shock location reasonably well, but the shock is 
weaker than that indicated by the experiment. In the aft region, the 
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lag entrainment solution agrees favorably with the data. Upstream, 
however, there is a discrepancy on both the upper and lower surfaces. 

The generally poor agreement between all the computed results and the 
experimental data is attributed to the limitations of small disturbance 

O 

theory. This conclusion is supported by a previous result for this 
same case where it was not possible to match the test data except by 
performing an inverse calculation. Thus, the good agreement obtained by 
the inviscid and wedge solutions on the upper surface ahead of the shock 
must be considered to be fortuitous. It is noted that the lift for 
these results must necessarily be in serious disagreement with the 
test. Corresponding viscous parameters for this case from the lag 
entrainment solution are shown in Figure 11. 

Using the steady-state solutions as initial conditions, the airfoil was 
oscillated sinusoidally in rotation about x = 0.4 with a reduced 
frequency k = 0.6 and an amplitude of 0.5°, duplicating one of the 
unsteady test conditions of Davis and Malcolm. The condition k = 0.6 
may violate the low frequency approximation, but the range of 
applicability of the governing equation can easily be extended by 

2 

including a high frequency term in the surface boundary condition . 

Three cycles of the pitching oscillation were computed with a value 
of At corresponding to 1° of rotation per time step. Instanteous 
pressure distributions for the last cycle of oscillation are shown in 
Figure 12. It is seen that the inviscid and wedge solutions vary only 
slightly from their respective steady-state values. The lag entrainment 
distribution, on the other hand, alters considerably during the course 
of the oscillatory cycle. This is due to the weakened shock which is 
present and allows for greater pressure fluctuations in the supersonic 
region. At t = 270 it can be seen that the pressure level crosses the 
sonic value (Cp*) more than once. Because this criteria is used to 
define and locate the wedge position of the viscous coupling, special 
provision had to be made to treat this case. This was done by forcing 
the code not to consider sonic locations upstream of a fixed value 
(x = 0.5). 
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A Fourier analysis of all solutions was performed on the upper surface 
pressure distributions for the last cycle of oscillation and the 
magnitude and phase of the first Fourier component was extracted. A 
comparison of these results is shown in Figure 13. In general the lack 
of agreement with the experiment is considerable. The inviscid solution 
fails to account for the pressure fluctuations ahead of the shock and 
also poorly predicts the shock location and pressure level in that 
region. The wedge solution has some of the same deficiencies, but 
compares better with the data near the shock and in the aft region. In 
the case of the lag entrainment result, the previously noted 
characteristics which were observed in the steady solution result in 
large fluctuations upstream of and at the shock. It appears that none 
of the computations have been able to closely reproduce the experimenal 
data. Once again this is attributed to the limitations of the theory. 

Corresponding instantaneous distributions of the boundary-layer 
parameters appear in Figure 14. Time histories of the instantaneous 
unsteady lift and moment coefficients and angle of attack are shown in 
Figure 15. The inviscid and wedge solutions are seen to be quite' 
similar. The lag entrainment solution differs from these not only in 
magnitude, but also in phase. This is not surprising based on the 
steady-state results of Figure 10. 

A second unsteady oscillating case was calculated for k = 0.4 and an 
amplitude of 1°. Pressure distributions are shown in Figure 16. The 
results of a Fourier analysis of the upper surface pressure distribution 
are compared with experiment in Figure 17. It is noted that the 
characteristics of the magnitude and phase are quite similar to the 
previous unsteady calculation. The inviscid and wedge solutions do not 
recover the upstream fluctuations and the inviscid result locates the 
shock downstream of its experimental value. While the lag entrainment 
result predicts upstream fluctuations, they are quantitatively 
inaccurate. In addiition, the phase in aft region does not agree with 
experiment. Corresponding instantaneous values of e and a* are shown in 
Figure 18 and time histories of the lift and momentum appear in Figure 
19. Phase differences between the respective solutions are apparent. 
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V. CONCLUSIONS AND DISCUSSIONS 


Two methods have been considered for including viscous effects in the 
compuation of unsteady transonic flows. One employs a simple phenomenological 
model as a computational artifice in accounting for the gross dominant effects 
of shock-boundary-layer interaction. The second, used in conjunction with the 
first, utilizes an additional set of integral equations to treat displacement 
effects of aft cambered regions and viscous trailing wakes witin the context 
of boundary-layer theory. Theoretical developments of the methods have been 
presented and modifications for their inclusion in the LTRAN2 code have been 
described. While the methods considered are relatively simple they are . 
thought to model the significant effects of viscosity and thus be useful for 
aeroelastic applications without degrading numerical efficiency. 

A number of computational results for both steady and unsteady flows have been 
generated. Comparisons with inviscid solutions and experimenal data have been 
made. Results which duplicate previous steady cases compare favorably, thus 
ensuring the integrity of modifications to the code. The unsteady 
calculations for the NLR 7301 airfoil exhibited serious deficiencies compared 
to experimental data. Although this result was not surprising, it is 
attributed to the severity of the case which may have violated some of the 
inherent theoretical assumptions. It is suggested that other unsteady 
computations be performed for a more thorough understanding of the limitations 
of the methods. 

All of the calculations presented here were carried out on the NASA AMES 
Research Center CDC 7600 computing system. Solutions for the wedge model were 
generated in virtually the same computational times as those for the inviscid 
LTRAN2 code. Inclusion of the lag entrainment equations required 
approximately 10 % additional computing time. Thus, the goal of maintaining 
efficiency appears to have been met. It is hoped that the modified code can 
thus extend the utility of LTRAN2 for the computation of unsteady transonic 
flows when viscous effects of shock-boundary-layer interaction, aft 
displacement, and trailing wakes are not negligible. 
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APPENDIX A 


At the leading edge of all blunt-nosed airfoils, the surface slope becomes 
infinite. Because this slope provides the flow tangency condition, the basic 
assumptions of small disturbance theory are violated in this region. Thus, 
since the solution in the immediate neighborhood of the leading edge cannot be 
described correctly by a numerical integration of the governing equation, a 
modification of the surface slope near the nose is often employed in transonic 
calculations. Such modifications are applied locally to enhance stability 
without degrading the results outside of the leading edge region. 


18 

One of the most common of these modifiations is the Riegels' rule . If 
i^- is the original airfoil surface slope, then the modified slope 
f^j^ is computed according to the following formula: 






This expression alters the slope in a smooth continuous manner at all points 
along the airfoil surface. When the slope is very large (i.e., near the 

leading edge) it is reduced to unity in the limit of f^ ♦ ». Elsewhere, 
when the slope is small, little change occurs. In the^Timit 

f* -► 0. f * = • Although the basis for Riegels' 

Aq AIiI Aw 

rule is incompressible theory, it has been widely used for computing steady 
transonic flows. For completeness, the modified and unmodified surface slopes 
for the NLR 7301 airfoil which were used for all the computations presented 
are given below. The unmodified slopes are those generated at computational 
mesh points by the LTRAN2 code which employs a cubic spline fit through 
airfoil coodinates for this purpose. 
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0.00250 

0.00750 

0.01250 

0.01750 

0.02250 

0.02750 

0.03250 

0.03813 

0.04563 

0.05500 

0.06500 

0.07500 

0.08563 

0.09875 

0.11703 

0.14204 

0.17188 

0.20313 

0.23438 

0.26563 

0.29688 

0.32813 

0.35938 

0.39063 

0.42188 

0.45313 

0.48438 

0.51563 

0.54688 

0.57813 

0.60938 

0.64063 


<0 

4 

X 1 

o 

^xm 

3.11794 

0.95222 

- 3.05200 

- 0.95029 

1.78953 

0.87295 

- 1.42208 

- 0.81800 

1.33984 

0.80140 

- 0.90111 

- 0.66942 

0.96358 

0.69387 

- 0.66624 

- 0.55445 

0.75954 

0.60485 

- 0.61575 

- 0.52432 

0.67158 

0.55752 

- 0.51744 

- 0.45956 

0.52192 

0.46269 

- 0.46085 

- 0.41854 

0.38352 

0.35809 

- 0.41565 

- 0.38381 

0.34466 

0.32585 

- 0.35792 

- 0.33699 

0.27566 

0.26575 

- 0.29837 

- 0.28591 

0.23412 

0.22796 

- 0.26036 

- 0.25196 

0.20572 

0.20150 

- 0.23001 

- 0.22416 

0.18233 

0.17937 

- 0.20005 

- 0.19617 

0.15831 

0.15637 

- 0.17457 

- 0.17197 

0.13440 

0.13321 

- 0.15341 

- 0.15164 

0.11412 

0.11339 

- 0.12556 

- 0.12458 

0.09122 

0.09084 

- 0.10286 

- 0.10232 

0.07817 

0.07793 

- 0.08716 

- 0.08683 

0.07083 

0.07065 

- 0.07088 

- 0.07070 

0.05150 

0.05143 

- 0.05092 

- 0.05085 

0.03627 

0.03624 

- 0.02748 

- 0.02747 

0.02288 

0.02287 

- 0.00449 

- 0.00449 

0.00028 

0.00028 

0.01371 

0.01370 

- 0.01312 

- 0.01312 

0.02928 

0.02927 

- 0.01191 

- 0.01191 

0.04426 

0.04422 

- 0.01294 

- 0.01294 

0.06119 

0.06108 

- 0.02988 

- 0.02987 

0.09409 

0.09368 

- 0.05300 

- 0.05293 

0.12610 

0.12511 

- 0.07613 

- 0.07591 

0.15442 

0.15262 

- 0.09763 

- 0.09717 

0.18383 

0.18080 

- 0.12661 

- 0.12561 

0.20768 

0.20334 

- 0.14333 

- 0.14188 

0.22189 

0.21662 
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X 


f'^ 

f" 

f" 


xo 

xm 

xo 

xm 


0.67188 

- 0.15362 

- 0.15184 

0.21845 

0.21342 

0.70313 

- 0.16744 

- 0.16514 

0.21495 

0.21015 

0.73438 

- 0.19029 

- 0.18694 

0.21649 

0.21159 

0.76563 

- 0.20953 

- 0.20508 

0.20620 

0.20195 

0.79688 

- 0.21820 

- 0.21319 

0.18921 

0.18591 

0.82813 

- 0.23116 

- 0.22522 

0.15643 

0.15455 

0.85797 

- 0.23315 

- 0.22706 

0.12979 

0.12871 

0.88235 

- 0.23304 

- 0.22696 

0.09068 

0.09031 

0.90000 

- 0.23274 

- 0.22668 

0.06194 

0.06182 

0.91500 

- 0.23253 

- 0.22649 

0.04077 

0.04074 

0.93000 

- 0.23214 

- 0.22613 

0.01876 

0.01876 

0.94500 

- 0.22968 

- 0.22385 

- 0.00639 

- 0.00639 

0.96000 

- 0.22468 

- 0.21921 

- 0.03459 

- 0.03457 

0.97500 

- 0.22097 

- 0.21576 

- 0.05940 

- 0.05930 

0.98875 

- 0.22569 

- 0.22015 

- 0.07582 

- 0.07560 

1.00000 

- 0.22791 

- 0.22221 

- 0.08597 

- 0.08566 
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APPENDIX B 


A brief description for use of the modified code is presented here. It is 
assumed that the user is familiar with the LTRAN2 program so that only the 
alterations governing its use will be considered. All user defined input 
parameters have assigned default values. The code will operate with any 
number of these redefined by the user as input. In many cases, if default 
values of some of the parameters are used, the resulting solutions may not 
have any practical significance. The following list itemizes the user defined 
parameters: 


Variable Name Variable Type Default Value Description 


Vise 

Logical 

FALSE 

Defines boundary-layer option 
FALSE— no viscous calculation 
TRUE— perform viscous calculation 

XOFFST 

Real 

0.02 

Offset of wedge from sonic point, 
^0* 

XPREC 

Real 

0.02 

Length of wedge precursor, Xp. 

XRAMP 

Real 

0.10 

Length of ramp, 

IBLCAL 

Integer 

1 

Iteration increment for 
performing viscous calculation 

GREEN 

Logical 

FALSE 

Defines lag entrainment option 
FALSE— use viscous wedge only 
TRUE — integrate lag entrainment 
equations 

REYINF 

Real 

1.0E7 

Chord Reynolds number based on 
freestrem conditions. Re 

00 



TINF 

Real 

300.0 

Freestream temperature in degrees 
Kelvin, T 

PRT 

Real 

0.9 

Turbulent Prandtl number, Pr^ 
(default value is for air) 

SO 

Real 

110.0 

Sutherland law viscosity constant 
in degrees Kelvin, So (default 
value is for air) 

IBLPRT 

Integer 

50 

Iteration increment for printing 
viscous parameters 


The logical variable VISC has been added to the namelist INPUT. If no viscous 
calculations are desired, VISC may be input as FALSE or the default value 
taken. This will permit the code to operate in the inviscid mode. If VISC is 
input as TRUE, it indicates that a viscous computation is to be performed. In 
addition the program expects to find another namelist of input variables 
immediately following INPUT. The new namelist name is VCALC. Note that even 
if no parameters are defined in VCALC, a card of the form $VCALC $ must appear 
in the input data. Conversely, if VISC is FALSE, the namelist VCALC must not 
appear. 

The remaining parameters are defined in the VCALC namelist. The variables 
XOFFST, XPREC, and XRAMP are the parameters defining the ramp geometry. In 
general the optimal choice for these is case dependent, but the default values 
have proven adequate for a large number, of both steady and unsteady 
computations. IBLCAL is an integer variable available to control how often 
the viscous displacement is updated. It is suggested that this be done 
every time step, but for certain applications an alternate procedure may be 
desirable. GREEN is a logical variable controlling selection of the lag 
entrainment calculation. If GREEN is FALSE only the viscous wedge is used. 
With GREEN set to TRUE the lag entrainment equations will be integrated. 
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When the lag entrainment option is selected, the variables REYINF, TIMF, PRT, 
and SO must be specified. These are the parameters which define the phj'sical 
flow under consideration. The default values will have no significance in 
this case. Finally, IBLPRT is an integer variable which determines the 
number of time steps between printing values of the viscous parameters. This 
printing gives distributions of 


1 


"n, Cg. and 

X 


5 * 

ir 


along the airfoil and trailing wake for the upper and lower surfaces.. Note 
that no printout occurs for the wedge solutions unless the lag entrainment 
equations are integrated. 

The input for running the code in its various modes is as follows. 

1. Inviscid The code will run inviscidly by setting 

Vise = .FALSE., or by using its default value. 

2. Wedge Only In namelist INPUT set VISC = .TRUE. Namelist 

VCALC must now appear even if default values are 
taken. XOFFST, XPREC, XRAMP, AND IBLCAL may be 
defined or default values taken. All other 
variables are not needed, but GREEN = .FALSE, if 
defined. This is also its default value. 

3. Lag Entrainment Set VISC = .TRUE, in namelist INPUT. In namelist 

VCALC, set GREEN=.TRUE. Ramp parmeters may be 
defined or default values taken but REYINF, TINF, 
PRT, and SO must be defined for the specific flow 
conditions. Other parameters are user optional. 
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In order to generate a steady viscous solution, the airfoil surface boundary 
condition must be augmented by the viscous displacement through an unsteady 
time integration. This is done by setting the amplitude of the unsteady 
motion (ALFONE) to zero. Choice of the time step is controlled by input 
variables XK and NSPC, which will be case dependent. Implementation of the 
wedge requies that the minimum pressure on the airfoil occur somewhere between 
the leading and trailing edge. If this situation is characterized by a shock, 
the wedge will be located within the shock profile. In addition, if the lag 
entrainment equations are to be solved, they will apply aft of the shock up to 
the downstream boundary. When no shock occurs, the wedge is not applied, but 
the lag entrainment solution (if this option is selected) will commence at the 
point of minimum pressure. 

Thus, to obtain the viscous solution the calculation should begin with a well 
behaved inviscid result as the initial condition. This result need not be 
converged, but should have the essential features of a converged solution. If 
it is not possible to obtain a reasonable inviscid result for the given flow 
conditions, these may be altered until the viscous modifications are 
implemented. For example, if in the inviscid solution the shock lies at the 
trailing edge, the Mach number and/or angle of attack may be reduced. This 
allows the shock to move upstream. Once the inviscid starting solution is 
obtained, viscous input parameters may be set to their desired values for the 
computation of the viscous solution. After the steady viscous solution has 
been generated, it will serve as the initial condition for any unsteady 
viscous calculation. 
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